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CNJ . We adapt the Smoothed-Particle Hydrodynamics (SPH) technique to allow a mul- 
• tiphasc fluid in which SPH particles of widely differing density may be freely infer- 
os | mixed. Applications include modelling of galaxy formation and cooling flows. 
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1 INTRODUCTION 

Since its introduction more than two decades ago by Lucy 
(1977) and Gingold & Monaghan (1977), Smoothed-Particle 
Hydrodynamics (SPH) has become one of the standard tech- 
niques for modelling astrophysical fluid flow (e.g. Evrard 
1988; Hernquist & Katz 1989, hereafter HK89; Thomas 
& Couchman 1992; Steinmetz & Muller 1993; Couch- 
man, Thomas & Pearce 1995, hereafter CTP95; Shapiro et 
al. 1996). SPH is fully Lagrangian, with the particles them- 
selves being the framework on which the fluid equations are 
solved, and so there is no grid to constrain the dynamic 
range or geometry of the system being modelled. This is of 
particular importance for phenomena involving the growth 
of gravitational fluctuations, such as cosmological structure 
formation (Bertschinger 1998) and star formation (Bhat- 
tal et al. 1998), and in an adaptive form (Wood 1981; Nel- 
son & Papaloizou 1994) the SPH algorithm will follow the 
wide range of densities encountered without difficulty. In 
contrast, computational demands severely restrict the dy- 
namic range of Eulerian finite-difference methods (e.g. the 
Piecewise-Parabolic Method of Collela & Woodward 1984), 
and to date the best three-dimensional Eulerian simulations 
of galaxy formation have a gas resolution of ~ 300 — 500 kpc 
(Blanton et al. 1999), although adaptive mesh refinement 
(AMR) techniques promise to greatly enhance the fixed- 
grid approach. SPH can be easily integrated with a range 
of A-body gravity solvers such as the tree algorithm used 
by HK89 and the Adaptive Particle-Particle, Particle-Mesh 
(AP 3 M) algorithm of Couchman (1991). 

However, SPH is not without its problems. The need 
for an artificial viscosity means that SPH resolves shocks 
poorly in comparison with finite-difference methods. Pair- 
wise artificial viscosities (Monaghan & Gingold 1983) gen- 
erally give the sharpest resolution of shocks, but also intro- 
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duce a large shear viscosity unless correction terms are used 
(Balsara 1995). There is some concern that these terms fur- 
ther degrade the shock capturing ability of SPH (Navarro & 
Steinmetz 1997; Thacker et al. 1998, hereafter T98). While it 
is possible to add a physical viscosity to SPH and solve the 
Navier-Stokes equations directly (e.g. Flebbe et al. 1994), 
the relatively simple SPH interpolation method is quite sen- 
sitive to particle disorder and tends to give large errors in 
the higher-order dissipative terms. Boundary conditions are 
also difficult to implement in SPH and do not sit naturally 
with the method. Generally boundaries are either periodic 
or situated far from the region of interest, and SPH is unsuit- 
able for simulations in which complex boundary conditions 
are of critical importance. 

Standard implementations of SPH also have a limited 
ability to resolve steep density gradients, and a number of 
numerical problems can occur when particles are close to, 
but not physically part of, a region of higher density. These 
arise because the usual formulation of SPH assumes that the 
density gradient across the smoothing kernel of each parti- 
cle is small. However, this is not true in many situations in 
which SPH is commonly used. In simulations of galaxy for- 
mation, for example, thermal instability causes cold, dense 
clumps to form within halos of hot gas, and density con- 
trasts of several orders of magnitude can occur within the 
typical smoothing length of the halo gas. Current implemen- 
tations of SPH will overestimate the density of the halo gas, 
leading to the gas cooling excessively and accreting on to 
the cold clump (T98; Pearce et al. 1999). Tittley, Couchman 
and Pearce (1999, hereafter TCP99) have also shown that 
the drag exerted on these protogalaxies by the intracluster 
medium can be seriously overestimated. These problems are 
thought to contribute to the overmerging commonly seen in 
simulations of structure formation including gas (Frenk et 
al. 1996). 

We wish to extend the method to cope with these prob- 
lems, along with multiphase fluids, such as the intraclus- 
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ter medium and cooling flows, which consist of an emulsion 
of discrete phases in which there is no correlation between 
the density of neighbouring particles. We should emphasise 
however that we are not seeking to model homogeneously 
mixed materials with different equations of state, such as 
dust and aerosols (Monaghan 1997). We make two changes 
to the standard implementation of SPH. Firstly, particles no 
longer have to try and maintain a fixed number of neigh- 
bours, as is required by the algorithm of HK89, and in- 
stead the smoothing length is adjusted to keep a density- 
weighted quantity constant (see Section 2/2). This prevents 
the smoothing length of low density particles from decreas- 
ing as they approach a high density region. Secondly we 
make the assumption that pressure, not density, is constant 
across the smoothing kernel. We then summate the local 
pressure at each particle and calculate the local density from 
the equation of state. This is a reasonable approach to take 
when the cooling time of the gas is greater than the local 
dynamic time, where a local pressure equilibrium can be ex- 
pected even when the local density gradient is steep. It will 
not be true in the presence of shocks (although neither is the 
assumption that density is smoothly varying) . As we will see 
and Bj, however, there is no evidence that 



in sections 1.1 



the shock capturing ability of SPH is degraded. 

The layout of the paper is as follows. In section y we de- 
tail our new implementation of SPH. In Section^ we present 
the results of a number of tests of the new method, and com- 
pare the performance with results produced using the SPH 
implementations of CTP95 and T98. This is followed in Sec- 
tion ^| by a discussion of our results. 



2 METHODOLOGY 

In SPH the fluid is represented by particles of known mass, 
mi, and specific energy, e^. Other properties must be inferred 
by averaging over a smoothing sphere that typically extends 
to enclose a fixed number, -/Vsph, of particles. For a contin- 
uous distribution, the average value of some quantity A at 
the location of particle i would be 



A(r) W(r-n)&V 



(1) 



where W is the smoothing kernel. In SPH the integral is 
replaced by a sum: 



3 



(2) 



that extends over all particles, j, within the smoothing 
sphere, pj is the density of particle j and so m,j/pj is its 
volume. Here 

W« = P"W(r«/ft«) I (3) 

where rtj is the separation of particles i andj, r%j = — rj|, 
and hij is the smoothing length. There are several possible 
forms for W . Throughout this paper, we use the standard 
form described in Thomas & Couchman (1992). The choice 
of hij determines the size of the region over which the density 
is to be averaged. Many authors take hij to be a symmetric 
function (for example the harmonic average) of h, and hj 
(see T99 and references therein). We prefer to set hij = hi. 



This has the advantages that we know in advance exactly 
how far we have to search and always have the same number 
of neighbours within the smoothing region. The integral of 
the smoothing kernel over all space is unity, which translates 
to the condition 
m,- , 



E 



-w. 



(4) 



although this will only be true in a statistical sense. 

Equation ^ seems to require the value of pj, which is 
not known in advance. However, it is possible to circumvent 
this by formulating SPH such that A is always a multiple 
of density and a known quantity. Taking A = p gives the 
standard SPH estimate for the density in the neighbourhood 
of a particle 



(5) 



Where expressions arise involving derivatives, we use 
the divergence theorem to transfer the derivative onto the 
smoothing kernel as follows: 

(VAi) = J VAWdV = - J AVWdV (6) 
(V.Ai) = J V.AWdV = -J A.VWdV (7) 

2.1 Evaluation of the density 

In multiphase SPH, it is important to distinguish between 
the density of a particle, pi, and the mean density in the 
neighbourhood of a particle, pi, defined by Equation^. The 
pressure of the gas is expected to be a smooth quantity be- 
cause (excluding shocks) the sound-crossing time is shorter 
than the flow time across the smoothing sphere. The density 
of individual particles, pi, can therefore be determined from 
an estimate of the local pressure and the equation of state, 



P = 2/3/96 



(8) 



The specific energy e is related to the gas temperature T 
by e = 3kT/2fj,m,H , where k is Boltzmann's constant, tuh is 
the mass of the hydrogen atom and fi — 0.6 is the relative 
molecular mass. The SPH estimate of the pressure is 



(Pi) = 2/3 



3 



miHiWi 



and the density of a particle can therefore be written 
3 (Pi) _ Y., "' ' .I 1 /, 



(Pi) 



2c, 



(9) 



(10) 



Variations in e, can cause large variations in density between 
neighbouring particles, whereas Equation [|, being an aver- 
age, varies only slightly between neighbours. Note that a 
cold, high-density clump of particles will contribute a lot of 
terms to Equation [H], but each of these will be given a low 
weight due to the presence of the ej term. 

The consistency condition, Equation ^j, becomes 



\ ^ 2m,jtjWij 



3 



(11) 



which will be true if the pressure is slowly-varying, as we 
have assumed. 
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2.2 The smoothing length 



It is common practice to use the following identity 



The SPH smoothing length is usually defined in terms of the 
radius of a sphere, centred on the particle in question, that 
encloses a fixed number of neighbours. One can either search 
for such a radius on each time-step, a time-consuming pro- 
cess, or allow the number of neighbours to vary and accept 
an estimate based on the number of neighbours found on 
the previous step: 



hi <— hi 



Q+(l- 



«)( 



iVs 



N, 



(12) 



Ni = 



where Ni is the actual number of neighbours and TVsph is the 
desired number (HK89). a is a convergence parameter: for 
a uniform distribution, choosing a = 1 returns the correct 
value of hi on the next step. Where large density contrasts 
are present, however, a must be reduced to avoid overshoot- 
ing and convergence is much slower. We use a convergence 
parameter a — 0.4 and Nsph = 32. 

For the extreme density contrasts we envisage, T98 have 
shown that Ni can oscillate between values which are alter- 
nately too low and too high. They solve this problem both 
by giving lower weight to particles at the extreme edge of 
the smoothing sphere, and by the introduction of a more so- 
phisticated convergence algorithm. In this paper, we suggest 
a simpler approach that uses a neighbour count weighted as 
follows: 

^Pr+p3 

3 3 

where pi is the mean density at particle i and the sum ex- 
tends over all neighbours within the smoothing sphere. Note 
that: 

(i) For a uniform distribution in which pi = pj, then 
Wij = 1 and Ni is simply equal to the number of neigh- 
bours. 

(ii) Particles do not notice neighbours with a much 
greater than average density (v>ij ~ for pi <C Pj). This 
prevents the common instability whereby changing the ra- 
dius of the smoothing sphere so as to include or exclude 
a high-density clump can dramatically change the number 
count. 

(iii) Neighbours with a lower than average density count 
at most double (wij ~ 2 for pi 3> Pj). While this could lead 
to an isolated cold particle in a hot medium having fewer 
neighbours than is desirable, such a situation is unlikely to 
arise in practice and Wij could be limited to a maximum 
value of unity if desired. 



2.3 The equation of motion 

In the absence of artificial viscosity, heating and radiative 
cooling, the equation of motion for a parcel of gas is simply 

dv _ VP 
dt ~~ ~p~' 

In the same spirit as above, we require an estimate of VP 
that does not depend upon the number density of particles. 
This is 



(14) 



(VP*) -^\/,,«,V.U".,. 



(15) 



VP 

p 



= V 



p 



+ 



p- 



Vp 



(16) 



to symmetrise the pressure force. This expression is not suit- 
able here, both because we are envisaging abrupt changes in 
density (so that Vp is undefined) and because the estimator 
for P/p 2 depends upon the density of the particles. Equa- 
tion lfij is obtained using the general identity 



i 



+ 



(17) 



p p" \P I p \p 

(Monaghan, 1992) and using a = 1 gives an alternative iden- 
tity that at first sight seems a little bizarre: 

VP VP P 
= + - VI. 

P P P 

This leads to the estimator for the force on particle i 



(18) 



^ 3 

3 



rairrij 



tjViWij tiViWji 



<p*> 



{Pi) 



(19) 



where we have chosen to use Wji in place of Wij in the 
second term on the right hand side in order to make the 
force symmetric. This gives 



<f*> = £ 



where 

2 ejViWij 
in — —rriim ■ 



(20) 



(21) 



3 <P*> ' 

The first set of terms in Equation ^o| is evaluated when cal- 
culating the SPH properties for particle i; the second set 
is accumulated in stages when calculating the properties of 
the neighbours. Note that the only change from the usual 
SPH formalism is that e; and tj have exchanged places in 
Equation [H]. However, this change means that the force on 
particle i is dependent only upon the local pressure gradient 
and not upon the density of its neighbours. 



2.4 Conservation of energy 

Conservation of energy is ensured by equating the heating 
to the PdV work done: 

J = --V.v. (22) 
dt p v ' 

By making use of the identity 

V.v = V.Av, (23) 



where Av = v — Vj, this can be put into a variety of forms. 
For example, 



de _ V.(PAv) Av.VP 
dt p p 



gives 

Ac, 
At 



2 \ ^ mjejVij.ViWij _ ... , 



3 ^ 

3 



(Pi) 



(24) 



(25) 



where Vy = Vj — Vj -the second term has vanished because 
it contains which is identically zero. Using 
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de _ PV.Av 

dt ~ p ' 

we obtain an alternative expression 



~dt) ~ 3 2.^ 



mjeiVij.ViW, 



dt/ 3^ (p~T 



— = / f ji- V J 



Combining Equations and ^7| we find that 



dc, 
~dt 



(26) 



(27) 



(28) 



where fy is the pairwise force given in Equation ^| 

Any of Equations ^7| or ^ (or any similarly derived 
equation) may be used as an estimator of de/dt (and hence, 
by rearrangement of Equation [22| for V.v). Indeed, if density 
is eliminated by substituting Equation |^ it is clear that the 
equations are identical, provided that our assumption that 
Pi and Pj are approximately equal is valid. Overall energy 
conservation is ensured as all three expressions give 

dE / de, 



dt 



dt 



+ Vi.fi 



(29) 



Under most conditions the performance of the three equa- 
tions is indistinguishable, and standard tests of SPH cannot 
tell them apart. However, when large temperature variations 
are present the double-sided form can lead to a large scatter 
in particle entropy. We use Equation ^ here. 

2.5 Artificial viscosity 

In the presence of shocks, we require a mechanism to con- 
vert relative motion into heat. In SPH, this is achieved via 
an artificial pressure, Q, that is added to the usual one in 
regions of convergent flow. This has the effect of replacing 
P by P + Q in Equations 14 an d [22 . In the current method, 
we replace fy in EquationsGfj and 



gij = f,,M,j[ -.VI, + a) 



with fij + gij , where 
(30) 



(c.f. Monaghan 
number 

0, 

Ma = { h ^j- 



Gingold 1983). Here the pairwise Mach 



oy (r 2 .+0.01h 2 ) ' 



> 
< 



the particle separation, n 
speed 



= r,- — n 



(31) 



and the sound 



Cj + c, 



where 

c= v/(10e/9). 



(32) 



(33) 



We adopt the values a = 1, /3 = 2. The use of cy in Equa- 
tion |3l] rather than Ci helps to limit the degree to which cold 
particles shock against dense clumps. 

In addition, a shear-correcting term (Balsara 1995) can 
be applied to limit the damping of shear flows. Following 
Steinmetz (1996) we use 



Mij -* Mij 
where 



| (V.v) 



(V.v)i| + |(V xv)i| + 0.Q0Qlci/hi 



(34) 



(35) 



For compressional flows k — 1 and Mi is unchanged, while 
in shearing flows k — > and the artificial viscosity vanishes. 
Only t he c osmological galaxy formation test described in 
Section 3.7 uses this correction term. 



2.6 Cooling 

In the presence of artificial viscosity and cooling, the Equa- 
tion of Energy Conservation becomes 



de 
dt 



P + Q 



V.v -5, 



(36) 



where £ is the emissivity (the emission rate per unit volume). 
The cooling function is interpolated from Sutherland & Do- 
pita (1993). To minimise problems of cooling during shock- 
heating (see Hutchings & Thomas 2000) , we allow the gas to 
cool only at the end of a time-step after the artificial viscos- 
ity term has been applied. The cooling is assumed to occur 
at constant density (the time-step ensures that this condi- 
tion is approximately satisfied), as described in Thomas & 
Couchman (1992). 



3 TESTS 

The tests discussed in the following sections were carried out 
using HydraQ (Couchman, Thomas and Pearce, 1995), an 
AP 3 M+SPH code modified to use the SPH method detailed 
in this paper. The simulations were performed on single- 
processor Sun and Intel workstations. 



3.1 Shock tube 

A standard test of gas dynamical codes is the Sod shock (Sod 
1978), which has been applied to SPH by many authors (e.g. 
Monaghan & Gingold 1983; Rasio & Shapiro 1991; T98). 
Analytic results are given by Hawley, Smarr & Wilson (1984) 
and Rasio & Shapiro (1991). This test is often carried out 
in one dimension, but this does not properly test particle 
interpenetration; we perform a three dimensional test, with 
the normal cubical simulation volume altered to match the 
geometry of a shock tube. Dimensions are 6 x 6 x 120 in code 
units and all boundaries are periodic. 

In a Sod shock two regions of gas with different den- 
sities are brought into contact, resulting in a shock wave 
propagating into the low-density gas and a rarefaction wave 
propagating into the high density gas. Between these two 
regions is a contact discontinuity, where the pressure is con- 
stant but the density jumps. Following T98 we use the initial 
conditions 



Pi - 

Pr 



4 

= 1 



P 

Pr 



1 

: 0.1795 



in 

V r 



0. for x < 
= 0. for x > 



(37) 



giving a shock Mach number M ~ 1.4. Both regions are 
allowed to evolve at constant temperature before being 
brought into contact, to allow the gas to relax to a physically 



t This code is in the public domain and can be downloaded from 
http: / /hydra, monaster. ca/| 
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Figure 1. Results from the Sod shock test after 30 time-steps. The variation in density, pressure, velocity and the entropic function 
A(s) = e/p( 7— 1 ' across the shock are shown. Analytic solutions are shown as solid lines, and the initial pressure and density profiles 
defined in equation Bjl are shown as dashed lines. The smoothing length in the vicinity of the shock is ~ 0.9 in code units. 



realistic initial state. A total of 7343 equal mass particles are 
used. 

Thirty timesteps are required to reach the state shown 
in Figure hi, by which time the shock front has moved around 
13 code units to the right. Both the shock and the contact 
discontinuity are broadened over a range Aa; « 3h. The 
results are in good agreement with the analytic solutions, 
although in common with other implementations of SPH 
the gradient of the rarefaction wave is too shallow. Flow is 
reasonably smooth, and there is no post-shock ringing. Fig- 
ure ^| compares the results from our code with those of T98 
and CTP95. There is very little difference between the three 
implementations, although the CTP95 implementation pro- 
duces a broader shock than the other two codes. This is due 
to the choice of artificial viscosity - CTP95 uses an artificial 
viscosity based on the divergence of the local velocity field, 
which is shown in T98 to be worse at shock capturing than 
the pairwise artificial viscosity used in T98 and our code. All 
three implementations model the rarefaction wave similarly, 
as viscous terms do not apply in expanding regions. 

With M^IA, this test represents a fairly weak shock. 
A strong adiabatic shock presents a more demanding test, 



as the pressure jump across the shock, and hence the Mach 
number, are infinite. In this case the jump conditions are 

pi = l Pi=0 v t = 1.333 for x < , . 

p r = 4 P r = 1.333 v r = 0.333 for x > 1 j 

although in our code we require particles to have a small 
minimum temperature to prevent numerical divergences in 
Equation [l(j, leading to Pi being slightly greater than zero 
and the Mach number remaining large but finite. The den- 
sity profile across the shock is shown in Figure ^. The new 
code is clearly capable of handling such shocks, and although 
some post-shock oscillation is visible it is not noticeably dif- 
ferent from the results obtained from the other codes. As the 
performance of our code is so close to that of T98 there is no 
evidence from shock-tube experiments that our assumption 
that pressure is smoothly varying has degraded the shock- 
capturing ability of the code. 



3.2 Adiabatic collapse 

One of the primary requirements for any hydrodynamics 
code that includes gravity is the ability to correctly follow 
the shock-heating of cold gas during gravitational collapse. 
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Figure 2. A comparison between the results from this code 
(dashed line) and those of CTP95 (dotted line) and T98 (dotted- 
dashed line). The shock front is at x ss 13. The broader shock 
produced by the CTP95 code is a result of the poor shock cap- 
turing of the artificial viscosity used. 
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Figure 4. Time evolution of the energies in the adiabatic col- 
lapse test. Results from the multiphase (solid line) and standard 
(dashed line) methods are shown, and the zero point is marked 
with a dotted line. The curves are, from top to bottom, the ther- 
mal, kinetic, total and potential energies. Normalised units are 
used. 




position 

Figure 3. The density jump across a strong adiabatic shock. 
Results from this code are shown. The smoothing length is ~ 0.8 
on the low-density side of the shock, and ~ 0.5 on the high density 
side, both in code units. 



A common test problem for SPH codes is the collapse of 
an initially isothermal sphere of gas (Evrard 1988; HK89; 
Steinmetz & Miiller 1993; T98), with an initial density pro- 
file 



p(r) = 



M(R) 1 



(39) 



2-kR 2 r 

To create this profile we translate particles radially from 
a uniform grid, which gives a lower sampling error than a 
random distribution of particles. Initial particle tempera- 
tures are set to the code minimum. A total of 8184 particles 
are used, with the gravitational softening set as 0.02R. Fol- 
lowing Evrard (1988) results are presented in normalised 
units, with density, internal energy, velocity, pressure and 



time normalised by 3M/4nR 3 , GM/R, (GM/R) 1/2 , pu and 
{R 3 /GM) 1/2 respectively. 

The gas in the sphere initially has negligible thermal 
energy, and collapses due to the lack of pressure support. 
During the collapse a central bounce occurs, causing a shock 
wave to propagate outwards through the gas, and the sphere 
should eventually reach virial equilibrium, with the ratio 
of thermal and gravitational energy approaching a value of 
—0.5. Figure |i| compares the evolution of the thermal, ki- 
netic, total and gravitational potential energies during the 
collapse for the multiphase and standard methods. The per- 
formance of the two codes is very similar throughout the 
bounce and subsequent expansion, with only relatively mi- 
nor differences apparent after time t ~ 1. 

Profiles of the system at time t = 0.8 are shown in Fig- 
ure H at which time the shock is located at r/r* w 0.2. 
The temperature jump across the shock is well modelled 
by both codes, and the post-shock conditions are very simi- 
lar. However, the density profile across the shock produced 
by the multiphase code is noticeably shallower. This is a 
result of the overestimation of the local pressure due to 
the steep pressure gradient in the shock, and is analogous 
to the smoothing of the density profile in the presence of 
steep density gradients that affects the standard implemen- 
tation of SPH. This error is unimportant in the force calcu- 
lations, which are dominated by the artificial viscosity dur- 
ing shocks, but may be significant when radiative cooling is 
implemented, as it leads to to excess cooling in the shock. 

In such circumstances it may be preferable to use the 
standard SPH estimate of the density pi for calculating the 
emissivity in regions where the local pressure gradient is 
steep. This switch can be implemented easily using either a 
pairwise condition based on the pressure of particles i and 
j or an SPH estimation of the local pressure gradient. Both 
methods effectively restrict the use of pi to the vicinity of 
the shock. When this approach is used, the density profile 



© 0000 RAS, MNRAS 000, 000-000 



Multiphase SPH 7 




across the shock in Figure g closely resembles that of the 
standard code, as would be expected, with the other profiles 
remaining unchanged. 



3.3 Density estimation in a two-phase medium 

In the standard formulation of SPH arbitrarily steep density 
gradients are smoothed over a distance representative of the 
local value of h. This is a cause for concern in simulations of 
cosmological structure formation, where cold dense clumps 
of gas - galaxies - form within diffuse halos of hotter gas and 
the density of halo gas can be overestimated by an order of 
magnitude or more in the steep density gradients around the 
largest galaxies. 

The smoothing of the density is a result of two ef- 
fects. Firstly, when Equation |^ is used, p tx h~ 3 for a fixed 
number of neighbours, and the estimated density becomes 
closely tied to the choice of smoothing length. Secondly, 
most adaptive forms of SPH update the smoothing length 
each timestep to try and keep the number of neighbouring 
particles approximately constant (e.g. the method of HK89). 
Particles act as tracers of the underlying mass distribution, 
and a high-density region will therefore contain more par- 



ticles than a region of lower density. When the HK89 algo- 
rithm is used, a particle close to a dense clump of gas will 
need to search little further than the edge of the dense region 
to find its required number of neighbours. If Equation |E] and 
the HK89 algorithm are used together then the estimate of 
the density of the particle becomes strongly dependent on 
its separation from the clump. 

Simply estimating the density using Equation [u] is not 
enough to cure this problem. While the high density par- 
ticles will be at low temperature in a region of pressure 
equilibrium, and will therefore contribute little to the esti- 
mate of the pressure, the size of the smoothing sphere, and 
hence weight given to neighbouring particles by the smooth- 
ing kernel Wij, still depends on the distance to the clump. 
It is therefore necessary to weight the neighbour count us- 
ing Equation |l3| so that a low-density particle close to a 
dense clump will assign a very low weight to particles in the 
clump and will search for neighbours as if the clump was not 
present. 

The benefits of our approach can be seen in Figure [j| 
Here we plot the change in smoothing length and the SPH 
estimate of the density as a hot gas particle moves past a 
cold, dense clump, grazing the surface of the clump at clos- 
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Figure 6. Change in the SPH estimate of the density p and the smoothing length h as a particle in a stream of hot gas moves past 
a cold, dense clump, grazing the surface of the clump at closest approach. The central density of the clump is 200 times that of the 
surrounding gas, with the temperature set so that the clump is in pressure equilibrium with its surroundings. Panel (a) shows results 
from the standard SPH approach, with density estimated using Equation |jj and the smoothing length updated using the method of 
HK89. These results are plotted using dashed lines in subsequent panels. Panel (b) shows the effect of our new /i-update algorithm 
when combined with the standard estimation of the density. Panel (c) shows the combination of the new estimate of the density given 
by Equation and the HK89 h-update method. Finally, panel (d) shows results from combining the new estimate of the density with 
the new h-update algorithm. 



est approach. The clump contains 420 particles, and has a 
central density 200 times that of the surrounding gas, with 
temperatures set so that the two phases are in pressure equi- 
librium. For the purposes of this test all inter-particle forces 
have been turned off, so that the particles move at constant 
velocity. Initially, the density of the hot particle is ~ 1.15 
and the smoothing length h ~ 0.75, both in code units. 

Panel (a) shows the results from the standard imple- 
mentation of SPH, in which the density is estimated using 
Equation |B| and the smoothing length is updated using the 
algorithm of HK89. Once the smoothing kernel of the hot 
particle overlaps the cold clump it finds many more than 
the desired number of neighbours and the HK89 algorithm 
starts to decrease the smoothing length in response, with h 
reaching a minimum of ~ 0.19 shortly after closest approach 
(this delay is due to the convergence parameter a in Equa- 
tion [l^ limiting the rate at which h can change). There is a 
corresponding increase in the estimate of the density, which 
peaks slightly earlier at p ~ 107, an overestimate of some 



two orders of magnitude. These values of h and p are typi- 
cal of cold particles on the surface of the clump, indicating 
that despite the large difference in temperature between the 
hot and cold particles the standard implementation of SPH 
treats them identically. 

Panel (b) shows the result of changing to the /i-update 
method described in Section [T^, with the density still esti- 
mated using Equation |H[ The peak density has only dropped 
slightly to p ~ 80, while the smoothing length remains fairly 
constant until shortly before the closest approach, when the 
density of the hot gas particle has increased to a point where 
particles in the clump become significantly weighted, leading 
to a decrease smoothing length. 

Panel (c) shows the combination of the new estimate of 
the density, given by Equation and the HK89 /i-update 
algorithm. The peak density is reduced to p ~ 10, which 
occurs when the smoothing length reaches it's minimum and 
the cold particles in the clump are weighted most strongly. 

Finally, panel (d) shows the results when the new es- 
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Figure 7. SPH estimation of the density across a density contrast 
of pi/pr = 100, which occurs at r = 0.108. Results from our new 
algorithm are shown as a solid line, while the results from the 
algorithm of T98 are shown as a dashed line. All units are code 
units. 




Figure 8. Average radial velocity after 10 time-steps across a 
density contrast Ph/pi = 100, which occurs at r = 0.108. Results 
from our code are shown as a solid line, and results from the code 
of T98 are shown as a dashed line 



timate of the density is combined with the new smoothing 
length update algorithm. The smoothing length remains vir- 
tually constant throughout the transit, with h decreasing to 
~ 0.72 shortly after closest approach and p only increasing 
from 1.15 to 1.61. 

The new /i-update algorithm does imply some addi- 
tional computational expense when strong density gradients 
are present. In the example presented here the hot particle 
has in excess of 450 neighbours when the smoothing sphere 
overlaps the whole of the cold clump, far more than the 32 
required by the HK89 algorithm. The overall computational 
expense will not be nearly as severe as this might suggest, 
however, as this idealised test follows a single particle next 
to a dense clump, where the weighting of neighbouring par- 
ticles is the most extreme. For a distribution of particles, 
such as the cosmological structure formation discussed in 
Section 3.7, the overhead is typically on the order of ~ 10% 



per timestep. 



3.4 Force estimation in a two-phase medium 

Figure Q shows the SPH estimate of the density profile across 
the boundary between two regions of different density, and 
demonstrates the sharpness with which density contrasts 
can be resolved by our method. To create the high density 
region we first evolve a cubical volume of gas at constant 
temperature, to ensure a relaxed particle distribution. We 
then extract a spherical region and compress it radially to 
achieve the desired density. This is then inserted back into 
the cubical simulation volume and the particles are allowed 
to move at constant radius from the centre of the box for a 
few time-steps to ensure a fully stable initial state. Here the 
sphere has a radius r = 0.108 in code units, and is 100 times 
more dense than the surrounding gas, with the temperature 
again set so that the two regions are in pressure equilibrium. 
Away from the boundary both methods correctly estimate 



the density, as we would expect. However, as discussed in 
the previous section, the standard implementation of SPH 
clearly overestimates the density close to the dense region, 
whereas the new algorithm gives a density contrast that is 
much sharper. 

In addition to the problems in estimating the density, 
unphysical forces can also occur in the presence of steep 
density gradients. In implementations of SPH which try to 
keep the number of neighbours constant, a particle near a 
density gradient will find many, if not all, of its neighbours 
in the region of higher density. The pressure gradient at the 
particle will therefore be highly asymmetric, leading, from 
Equation |l4|, to a force that acts to push the particle away 
from the high density region. The strength of this force will 
depend on the magnitude of the density contrast, saturating 
once the particle finds all its neighbours in the high density 
region, and a dense clump of gas will therefore rapidly empty 
the surrounding region of particles. 

We can see the effect of this asymmetric pressure gra- 
dient in Figure ^| where we plot the radial velocities after 
10 time-steps of particles in or near the dense clump. The 
scatter in velocity due to thermal motion of the hot gas is 
around ±0.05 in both cases. The standard SPH code pro- 
duces a large outward velocity in the hot gas, which evac- 
uates a space between the two phases within a few tens 
of time-steps. This effect can be clearly seen in Figure 
where the number density of particles drops to zero between 
0.108 < r <0.13. While our method also produces excess ve- 
locities at the boundary, implying that the spurious pressure 
gradient has still not been completely eliminated, the mag- 
nitude of the outward velocity has been greatly reduced and 
the effect is far less systematic. 



3.5 Drag 

Drag resulting from gas dynamical forces is an important 
factor in simulations of cosmological structure formation, as 
incorrectly estimating the drag can bias both the distribu- 
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Figure 9. Number density of particles per unit volume after 25 
time-steps across a density contrast phi Pi = 100, which occurs at 
r = 0.108. Results from our code are shown as a solid line, and 
results from the code of T98 are shown as a dashed line. The line 
marking the boundary between the two regions has been omitted 
for clarity. The empty region around the dense central clump can 
be clearly seen. 
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Figure 10. Rate of acceleration of a dense clump in a flow with 
Mach number 0.5. Results from our code are shown as a solid 
line, results from T98 are shown as a dashed line, results from 
CTP95 are shown as a densely dotted line and the prediction of 
Equation ^ are shown as a sparsely dotted line. All codes produce 
an excess drag. 



tion of matter in clusters and the size and number of objects 
formed. Frenk et al. (1996) have suggested that excessive 
drag can worsen the overmerging problem seen in N-body 
simulations, and TCP99 have shown that standard imple- 
mentations of SPH can significantly overestimate the drag 
on a cold clump of gas moving through hot gas representa- 
tive of the intracluster medium. This excess drag is caused 
by accretion of gas from the ICM on to a shell around the 
clump, where it is held by forces arising from the miscalcula- 
tion of the pressure gradient around the clump (a discussion 
of this effect can be found in TCP99). This accretion leads 
to an increase in the effective radius of the clump, and hence 
the drag. This drag is most severe when the clump is mov- 
ing subsonically, and at higher velocities the accretion of gas 
decreases, largely vanishing when the Mach number M > 2. 

In order to examine the drag introduced by our imple- 
mentation of SPH, we consider the case of a clump of cold 
gas initially at rest in a stream of fast moving, diffuse gas 
(this is identical to the method of TCP99, except we work 
in the rest-frame of the cold clump) . Assuming collisions are 
inelastic, the clump will be accelerated to the flow velocity 
vo at a rate 



v = v 



+ 1 



where k is a constant 
2nr 2 p v 



k = 



M 



(40) 



(41) 



in which po is the density of the diffuse gas, vo the flow 
velocity, r is the radius of the clump and M is the mass of the 
clump (all in code units). For our tests M = 50, r = 0.5, the 
velocity and density of the flow are set to unity and the Mach 
number of the flow is determined by setting the temperature 
of the gas as required. We use a Mach number M ~ 0.5 here, 



which is well inside the regime in which TCP99 show that 
drag becomes excessive. Figure [H] compares the results from 
our new code with those obtained the codes of CTP95 and 
T98 (the code used by TCP99 is essentially the same as 
that used in T98). The results of our code are clearly an 
improvement, being close to the prediction of Equation 
This is due to the elimination of the accretion that occurs 
in the other codes, as the asymmetric pressure gradients 
causing accretion are greatly reduced by our method. The 
difference between CTP95 and T98 is mainly a result of the 
kernel smoothing used in T98, which is shown in TCP99 to 
further increase the effective cross section of a dense clump. 
At higher Mach numbers the differences between the codes 
are less pronounced as the accretion decreases once M > 1, 
although T98 continues to give the largest drag of the three 
codes. 



3.6 Multiphase cooling flows 

Clusters of galaxies contain large quantities of hot, X-ray 
emitting gas, often concentrated around the most massive 
galaxy in the cluster. In ~ 70% — 90% of cases, gas in the 
centre of the cluster has a radiative cooling time less than 
the Hubble time, Hq 1 (Edge, Stewart & Fabian 1992; White, 
Jones & Forman 1997), and as this gas cools surrounding gas 
will move inwards to maintain pressure support, initiating a 
large-scale motion known as a cooling flow (see Fabian 1994 
for a review). The mass deposition rate can be as high as 
10 3 Moyr -1 (Fabian et al. 1985), and is often observed to 
vary with radius roughly as M(< r) tx r within the radius 
'"cool in which the cooling time is less than a Hubble time 
(e.g. Thomas, Fabian and Nulsen 1987, hereafter TFN87). 
This is generally taken as evidence for a multiphase flow, in 
which thermal instability is causing the denser gas to cool 
out of the flow at larger radii. 

Modelling a multiphase flow is impossible with the usual 
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implementation of SPH as density is a locally averaged quan- 
tity, making the flow inherently single phase. In contrast, 
our method allows a wide range of densities, provided that 
a local pressure equilibrium exists. This is a reasonable as- 
sumption for particles with temperatures above 10 6 K, which 
will have cooling times much greater than the local sound 
crossing time (Nulsen, 1988). Particles below 10 6 K will cool 
to 10 4 K within a few time-steps, at which point they are 
removed from the flow. 

To test the ability of our code to model a fully 
multiphase environment, we examine a simple constant- 
pressure, spherically-symmetric cooling flow. The distribu- 
tion of phases in the flow is described by the fractional vol- 
ume distribution f(p,r) introduced by Nulsen (1986, here- 
after N86), where /dp is the fractional volume occupied by 
phases in the density range p to p+dp. N86 considered many 
analytic forms for /. Here, we take 

( 3 - a) / P_ 
Po I P° , 



f(fi,r) = 



P > Pa 



(42) 



with po being a minimum density and a the temperature de- 
pendence of the cooling function A tx T a . For the purposes 
of this test we replace the Sutherland & Dopita (1993) cool- 
ing function with a pure power law in which a = 0.5. Then 
Equation ^ can be integrated (N86) and gives a mass de- 
position profile 



M oc r v 

where 
3(2 



9 



i] 



(43) 



(44) 



(5 -2a) 8' 

Particles are placed randomly within a cubical simula- 
tion of volume (200kpc) 3 , and are then allowed to evolve 
at constant temperature until spurious fluctuations arising 
from the initial particle distribution have died away. Par- 
ticles are then ordered in terms of their distance from the 
centre of the box, and translated radially so as to match the 
mass profile given by 

M tx po(r)r 3 (45) 
where the density profile 

po(r) oc r~ T^^T oc r~ 3/4 (46) 

is derived in TFN87. Particles translated outside the bounds 
of the box are discarded. Particle densities are drawn at ran- 
dom from the volume fraction distribution given by Equa- 
tion with particle temperatures set so as to maintain con- 
stant pressure. Here, the outer temperature is T ~ 5 x 10 7 K , 
the inner temperature T ~ 10 6 _ftT, the central density is 
0.1cm -3 and the average outer density is ~ 5 x 10 _3 cm -3 . 
A total of 20252 particles are used. 

Figure [ll] shows the mass deposition profile M(< r), 
produced by the multiphase method after 500 timesteps. 
Roughly 11% of the gas has cooled from the flow. The 
line marks a least-squares fit to the data, with slope r\ = 
1.26±0.01. The points within 15kpc of the centre of the flow 
have been excluded from this fit, as there are too few parti- 
cles here for the mass deposition profile to be well sampled, 
as have particles with r > lOOkpc which lie in the corners of 
the cubical simulation volume. The slope is slightly steeper 
than the theoretical index of r\ = 9/8. Figure ha shows the 




10 

radius (kpc) 

Figure 11. Mass deposition profile produced by the multiphase 
code after 500 timesteps. Approximately 11% of the gas has 
cooled out of the flow. The solid line represents a least squares 
fit to the profile between 15kpc < r < lOOkpc, and has a slope 
1.26 ±0.01. 




P/<P> 

Figure 12. Mass distribution function pfdp of the uncooled gas 
at radii 50kpc < r < 60kpc for the multiphase code. We would ex- 
pect the distribution to be unchanged from the initial distribution 
given by Equation with the points on a power-law oc p~ 2 5 . 
The dashed line represents a least-squares fit (excluding the first 
two points), and has a slope —2.75 ± 0.06. 



mass distribution function pfdp for particles at radii of be- 
tween 50 and 60 kpc. In theory we would expect the distri- 
bution of phases given by Equation ^ to remain unchanged 
with time, giving pfdp oc p -2 ' 5 . The dashed line represents 
a least-squares fit to the data, with gradient —2.75 ± 0.06. 
This is again steeper than we would expect. This is proba- 
bly due to the phases not comoving, an assumption made in 
N86. There is no condition in our code to enforce this, and 
Equation [l4| implies that high-density particles will receive 
a smaller pressure force than low-density ones. The artificial 
viscosity will limit the degree to which particles can inter- 
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Figure 13. Mass deposition profile produced by the standard 
implementation of SPH after 500 timcstcps. Approximately 5% 
of the gas has cooled out of the flow. Virtually all the matter has 
been deposited within the central ~ 8 kpc, as would be expected 
of a single-phase cooling flow. 




P/<P> 

Figure 14. Mass distribution function p/dp of the uncooled gas 
at radii 50kpc < r < 60kpc for the standard implementation of 
SPH. All particles have densities close to the mean, and there is 
no trace of the original density distribution. 



penetrate, ensuring that the flow is largely comoving, but it 
may not be sufficient to ensure that no slippage occurs. Ap- 
plying a large bulk viscosity term in Equation |3^ forces the 
particles to comove, and flattens both slopes towards their 
theoretical values. However, this is not suitable for general 
application as it degrades the shock capturing ability of the 
method. 

The mass deposition profile produced by the standard 
implementation of SPH is shown in Figure [l3| 500 timesteps 
have passed and this time about 5% of the gas has cooled. 
In contrast with the multiphase results, this implementation 
produces a mass deposition profile that is centrally con- 
centrated, as would be expected from a single-phase flow. 



This is supported by the mass distribution function, which 
is plotted in Figure U4. The gas has clearly evolved back to 
a single-phase state, with all particles having a density close 
to the mean, and no trace of the original density distribu- 
tion remains. In addition, the gradient of the density profile 
has steepened from an initial value of p oc r -0 ' 75 (as given 



close to the theoretical 



by Equation ^(j) to p oc r~ 
single-phase gradient pec r~ 12 (Thomas, 1988). 

When compared to the standard implementation of 
SPH, it is clear that our method performs well. Both the 
mass deposition profile and the distribution of densities in 
the flow are close to the theoretical values, whereas the stan- 
dard implementation fails to reproduce the properties of the 
flow correctly. 



3.7 Cosmological galaxy formation 

Our final test involves a simulation of the formation of a 
cluster of galaxies. We use the initial conditions from the 
fiducial simulation of Kay et al. (2000), who used it to test 
in detail the effect of varying a number of numerical and 
physical parameters. The simulation uses the standard cold 
dark matter (SCDM) cosmology, with f2 = 1, A = and 

h = 0.50. The baryon fraction was set from primordial nu- 
cleosynthesis constraints, Q.bh 2 = 0.015 (Copi, Schramm & 
Turner 1995), and an unevolving gas metallicity of 0.5Zq 
was used. The initial fluctuation amplitude was set so that 
the model produces the same number density of rich clus- 
ters as is observed today, with as, the present-day linear 
rms fluctuation on a scale of 8/i _1 Mpc, set to 0.6 (Eke, 
Cole & Frenk 1996; Vianna & Liddle 1996). Until z ~ 1 
we adopt a comoving /3-spline gravitational softening equiv- 
alent to a Plummer softening of 20/i _1 kpc, after which it is 
switched to a fixed physical softening of ~ 10 /i _1 kpc. The 
minimum SPH resolution is set to match the constraints im- 
posed by the gravitational softening. 32 3 particles of both 
dark matter and gas were used, giving a mass per particle 
of 8 x lO 9 ft _1 M for the dark matter and 5 x 10 8 ft _1 M© 
for the gas. The simulation was started at z ~ 24 and was 
evolved to the present day. 

For the purposes of this simulation, Equation ^ incor- 
porates a source term 7i to reflect the heating of gas due to a 
photoionising background. We assume that the background 
has a standard power-law form 

J{v) = J 2 i(z) x 10~ 21 ( — ) ergss* 1 cm" 2 sr" 1 Hz" 1 ,(47) 



where 



Jai(«) = 



T° 
' >'? i 



(5/(1- 



(48) 



is the flux at the Hi Lyman limit (Navarro & Steinmetz, 
1997); we take J21 = 1 and a — 1 here. Photoheating is 
implemented following Theuns et al. (1997). The ultraviolet 
background has the effect of imposing a minimum tempera- 
ture on the gas, rapidly heating it to ~ 10 4 K and ensuring 
that pressure gradients in the gas remain shallow. If the 
effects of photoionisation are not included, then problems 



$ Ho = 100/ikms _1 Mpc" 
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Figure 15. The temperature-density phase diagram of the gas 
at z = in the multiphase cluster-formation simulation. Density 
is in units of the mean density. The diagram is divided into two 
phases; a galaxy phase, containing gas which is overdense by at 
least a factor of 500 and has a temperature of T < 10 K and a 
hot phase, containing all gas with a temperature of T > 10 5 K 
and the gas which is not sufficiently overdense to be considered 
part of the galaxy phase. 



can occur when particles which have been in free expan- 
sion since the start of the simulation, cooling to very low 
temperatures with T oc (1 + z) 2 , encounter the accretion 
shock at the outer edges of halos. This shock is generally 
poorly-resolved, and neighbouring particles that have passed 
through the shock can make an overwhelming contribution 
to the estimate of the pressure of the cold particle, which 
can lead to the density being overestimated, potentially by 
several orders of magnitude. Incorporating photoionisation 
serves as a simple way to limit this effect; an alternative is 
to use the standard SPH estimate of the density for calcu- 
lating the emissivity in regions whe re the pressure gradient 
is steep, as discussed in Section 3.2 

Figures [l5] and |l6| illustrate the temperature-density 
distribution of baryonic matter at z — produced by the 
multiphase and standard codes. We have divided the gas 
into two phases; a galaxy phase containing gas which is 
overdense by at least a factor of 500 and has a tempera- 
ture of 10 3 K< T < 10 5 K and a hot phase which includes 
all gas with a temperature of T > 10 K and the gas with 
T > 10 3 K which is not sufficiently overdense to be consid- 
ered part of the galaxy phase. Almost all the particles in the 
galaxy phase lie along a line with T ~ 1 — 2 x 10 4 K, which 
marks the point at which the radiative cooling and pho- 
toheating rates balance. Most of these particles are in the 
form of dense clumps with a size similar to the gravitational 
softening length - we refer to these clumps as galaxies. The 
properties of these galaxies are calculating by first extract- 
ing all particles within the galaxy phase from the simulation 
volume and then running a 'friends-of-friends' group finder 
(Davis et al. 1985). The most significant difference between 
the two simulations is in the mass of the largest galaxy, 
which is nearly 50% more massive in the single-phase simu- 
lation. Kay et al. (2000) find that this is a result of excessive 
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Figure 16. The temperature-density phase diagram of the gas at 
z = in the standard SPH cluster-formation simulation. Density 
is in units of the mean density. The phases are defined in the same 
way as in Figure Il5l 
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Figure 17. Radial temperature profile through the hot gas halo 
surrounding the largest galaxy found in the multiphase simula- 
tion. 



cooling, as decoupling the hot and cold phases (Pearce et 
al. 1999) reduces the final mass of the galaxy to a more rea- 
sonable value. 

The hot phase consists of particles that have been 
shock heated during gravitational collapse, with the bulk 
kinetic energy of the gas being converted to heat. Figure |l5| 
shows a significant quantity of hot (T > 10 7 K), high-density 
(p/(p) > W 4 ) gas which is not present in Figure |l6| but is 
seen in the single-phase simulations when radiative cooling 
is not allowed. This gas is located near the centre of large 
dark matter halos, close to the central galaxy, and overes- 
timation of the already high gas density by the standard 
method leads to the gas rapidly cooling and being accreted 
by the galaxy. In contrast, the multiphase method correctly 
estimates the density of the gas, resulting in a slower cooling 
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Figure 18. Radial temperature profile through the hot gas halo 
surrounding the largest galaxy found in the single phase simula- 
tion. 



rate and the gas remaining at high temperatures throughout 
the halo. This effect can be seen in more detail in Figures [l?] 
and jl8| which show the radial temperature profile of the gas 
around the most massive galaxy found in the simulation vol- 
ume. The size of the galaxy is on the order of the softening 
length, ~ 10ft _1 kpc, and is not shown here. The results from 
the standard implementation clearly show two of the prob- 
lems inherent in the method. Firstly, particle temperatures 
drop sharply within r < lOOkpc - this is clear evidence for 
overcooling. No such effect is visible in Figure [l7| and the 
lack of overcooling in our method is probably the principle 
reason for the difference in mass of the largest objects in the 
two simulations. Secondly, there is a complete absence of 
particles at radii r < 50kpc, whereas particles are found all 
the way in to r ~ lOkpc in the multiphase simulations. This 
is an example of the effect examined in Section 3.4, with 



particles close to the central galaxy being forced away by an 
artificially asymmetric pressure gradient in the single-phase 
method. 



4 DISCUSSION AND CONCLUSIONS 

We have presented a multiphase implementation of 
Smoothed-Particle Hydrodynamics (SPH), along with a 
number of tests to compare the performance of our method 
with standard implementations of SPH. The usual SPH for- 
malism assumes that density is a smooth quantity, varying 
negligibly on distances on the order of a typical smoothing 
length. This is clearly not true in many situations in which 
SPH is applied, such as simulations of galaxy formation, in 
which large density contrasts are present. However, the pres- 
sure of the gas is expected to be a much smoother quantity 
because in almost all situations the sound-crossing time is 
shorter than the flow time across the smoothing sphere. We 
therefore summate the local pressure at each particle, and 
calculate the density from the equation of state. 

One situation in which our assumption will definitely 
not be true is in the presence of shocks, in which density 



sure. Sections 3.1 and 3.2 demonstrate that our method 
handles shocks acceptably, although the estimate of the den- 
sity in the shock can be inaccurate when strong shocks are 
resolved poorly. This does not alter the force calculation, 
which is dominated by the artificial viscosity under such con- 
ditions, but is potentially significant when radiative cooling 
is implemented, and it may be may be preferable to use the 
standard estimate of the density for calculating the emissiv- 
ity in regions where the local pressure gradient is steep. In 
Section 3.3 we examine the degree to which steep density 
gradients can be resolved by the two methods. Standard im- 
plementations of SPH are shown to severely overestimate the 
density of particles close to a region of high density, while 
these particles have their density correctly est imated by our 
method. In addition, we show in Section 3.4 that unphys- 
ical forces can occur in the presence of steep density gra- 
dients, and while such forces are not completely eli mina ted 
by our method, they are greatly reduced. Section 3.5 ex- 



amines the drag introduced by our implementation of SPH. 
We find that at low Mach numbers drag is greatly reduced 
compared to the SPH implementations of CTP95 and T98, 
while at higher velocities the results are in broad agreement 
with the findings of TCP99, depending mainly on the choice 
of smoothing kernel symmetrization. 

In Section we examine a simple spherically- 

symmetric, constant pressure cooling flow. The fractional 
volume distribution /dp of phases within the flow is taken 
to be a pure power-law, which is shown by N86 to remain 
unchanged with time and deposit mass as M(< r) tx r 9//8 . 
The standard implementation of SPH is shown to be unable 
to reproduce the expected behaviour of this cooling flow, 
with conditions rapidly returning to a single-phase state. 
Our code performs well, although the gradient of both the 
mass distribution pfdp and the mass deposition profile M 
are too steep. This is probably due to phases not comoving in 
the flow, which is an assumption made in N86. Application 
of a large bulk viscosity to force particles to comove appears 
to reduce the problem, although this is not a suitable for ap- 
plication generally as it results in the shock capturing ability 
of the code being significantly degraded. 



In Section 3.7 we examine the formation of a cluste r of 
gala xies . The idealised problems examined in Sections 3.3 
and 3.4 are shown to have real analogues in the simulation 



using the standard SPH implementation, with overcooling 
resulting from overestimating the density of halo gas be- 
ing clearly visible, and halo gas being forced away from the 
central galaxy. No such effects are visible in the simulation 
using our method. The lack of overcooling is also apparent 
in the masses of the galaxies formed, with the largest galaxy 
being nearly 50% more massive in the single-phase simula- 
tion, in agreement with the findings of Pearce et al. (1999), 
who found that decoupling the galaxy from the hot halo gas 
produced a similar effect. 

Our method is an alternative to the standard formula- 
tion of SPH. In simulations without large density contrasts, 
the two give very similar results. However it represents a 
significant improvement over the standard implementation 
of SPH when the gas component cannot be assumed to be a 
single-phase fluid, such as galaxy formation and cluster for- 
mation. In addition, fully multiphase fluid flow can be mod- 
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elled, allowing SPH to be applied to simulations of cooling 
flows and the intracluster medium. 
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